function [M, pStar,exitflag] = get_pstar(V, Pgrid, ~)
% MAX value and optimal price at each state with quadratic interpolation on V

[nump]=length(V);
M = NaN; 
pStar = NaN;   

[~, maxind] = max(V);                                    
if  min(maxind)==1 || max(maxind)==nump                                     % make sure solution is interior
    disp ('Corner solution. Increase grid range')
    exitflag=0;
else
    localx = [Pgrid(maxind-1)';Pgrid(maxind)';Pgrid(maxind+1)'];         % 3 points on price grid around optimum
    XMATstack = [ones(size(localx));localx;localx.^2];                   % regressor matrix: const, grid and grid^2
    XMAT = reshape(XMATstack(:),3,3);                            % regressor matrix
    localV = V(maxind-1:maxind+1);                     % explained variable
    betacoeff = (XMAT'*XMAT)\XMAT'*localV;                           % regression coefficient
    pStar = -0.5*betacoeff(2)/betacoeff(3);                     % maximant of interpolant
    M  = [1 pStar pStar^2]*betacoeff;                 % maximum of interpolant
    exitflag=1;
end